{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ordinal Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some data are discrete but intrinsically **ordered**, these are called [**ordinal**](https://en.wikipedia.org/wiki/Ordinal_data) data. One example is the [likert scale](https://en.wikipedia.org/wiki/Likert_scale) for questionairs (\"this is an informative tutorial\": 1. strongly disagree / 2. disagree / 3. neither agree nor disagree / 4. agree / 5. strongly agree). Ordinal data is also ubiquitous in the medical world (e.g. the [Glasgow Coma Scale](https://en.wikipedia.org/wiki/Glasgow_Coma_Scale) for measuring neurological disfunctioning). \n",
    "\n",
    "This poses a challenge for statistical modeling as the data do not fit the most well known modelling approaches (e.g. linear regression). Modeling the data as [categorical](https://en.wikipedia.org/wiki/Categorical_distribution) is one possibility, but it disregards the inherent ordering in the data, and may be less statistically efficient. There are multiple appoaches for modeling ordered data. Here we will show how to use the OrderedLogistic distribution using cutpoints that are sampled from Improper priors, from a Normal distribution and induced via categories' probabilities from Dirichlet distribution. For a more in-depth discussion of Bayesian modeling of ordinal data, see e.g. [Michael Betancourt's Ordinal Regression case study](https://betanalpha.github.io/assets/case_studies/ordinal_regression.html) \n",
    "\n",
    "**References:**\n",
    " 1. Betancourt, M. (2019), “Ordinal Regression”, (https://betanalpha.github.io/assets/case_studies/ordinal_regression.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install -q numpyro@git+https://github.com/pyro-ppl/numpyro"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "from jax import numpy as np, random\n",
    "\n",
    "import numpyro\n",
    "from numpyro import handlers, sample\n",
    "from numpyro.distributions import (\n",
    "    Categorical,\n",
    "    Dirichlet,\n",
    "    ImproperUniform,\n",
    "    Normal,\n",
    "    OrderedLogistic,\n",
    "    TransformedDistribution,\n",
    "    constraints,\n",
    "    transforms,\n",
    ")\n",
    "from numpyro.infer import MCMC, NUTS\n",
    "from numpyro.infer.reparam import TransformReparam\n",
    "\n",
    "assert numpyro.__version__.startswith(\"0.19.0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Generation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, generate some data with ordinal structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "value counts of Y:\n",
      "1    19\n",
      "2    16\n",
      "0    15\n",
      "Name: Y, dtype: int64\n",
      "mean(X) for Y == 0: 0.042\n",
      "mean(X) for Y == 1: 0.832\n",
      "mean(X) for Y == 2: 1.448\n"
     ]
    }
   ],
   "source": [
    "simkeys = random.split(random.PRNGKey(1), 2)\n",
    "nsim = 50\n",
    "nclasses = 3\n",
    "Y = Categorical(logits=np.zeros(nclasses)).sample(simkeys[0], sample_shape=(nsim,))\n",
    "X = Normal().sample(simkeys[1], sample_shape=(nsim,))\n",
    "X += Y\n",
    "\n",
    "print(\"value counts of Y:\")\n",
    "df = pd.DataFrame({\"X\": X, \"Y\": Y})\n",
    "print(df.Y.value_counts())\n",
    "\n",
    "for i in range(nclasses):\n",
    "    print(f\"mean(X) for Y == {i}: {X[np.where(Y == i)].mean():.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABDi0lEQVR4nO3deXhU5f3//+d9ZiazZCWQhH3fRMqiEayiolYEUVvqrj+1FUWqLCpaBcSlKCAqq1utWn7afmot1mqlCrigiOwQdgKEkBAgIQvZZzLLub9/DMEFErLMzJnJ3I/r4pIkwzkvMzPnPec+93nfQkqJoiiKEn00owMoiqIoxlAFQFEUJUqpAqAoihKlVAFQFEWJUqoAKIqiRCmz0QEao02bNrJr165Gx1AURYkomzdvLpJSpvz8+xFVALp27cqmTZuMjqEoihJRhBA5Z/q+GgJSFEWJUqoAKIqiRClVABRFUaKUKgCKoihRShUARVGUKKUKgKIoSpRSBUBRFCVKGV4AhBAmIcRWIcSnRmdRlKbSdR3VWl2JNIYXAGAysMfoEIrSVJ9++imXX345V189gmPHjhkdR1EazNACIIToCIwG3jIyh6I0x549e5BS4nLVkJ2dbXQcpZE2b97M8uXL2bVrl9FRQs7oVhALgD8C8XU9QAgxDhgH0Llz59CkUpRGKC4uJtYiqfIIiouLjY6jNEJpaSmPPPIIUkpirDGsWL4CTQuHgZHQMOz/VAhxLXBcSrm5vsdJKd+UUqZLKdNTUk7rZaQohis8XkD3eC8CKCwsNDqO0gj5+flIKZEpEneNm5KSEqMjhZSRpe5i4HohxCHgfeAKIcTfDMyjKE1SUFBAil0nySYoKCgwOo7SCEePHgVAtpU/+TpaGFYApJRTpZQdpZRdgVuBr6SU/59ReRSlKaqrqymvqCTF7qONzcuxKDuARLq8vDzghwJw5MgRI+OEXPQMdilKENQeQNLsOmk2L0fyDhucSGmM3NxctFgNEgANDh06ZHSkkAqLAiClXCWlvNboHIrSWKcKgEMnzaFTWFyCy+UyOJXSUAeyDuCL84EGIkGoAqAoSsPl5OQggHYOH+1jfYD/U6US/jweD7k5ucgk//CPL8FH5r5Mg1OFlioAitIMOTk5pDggxgQdThaAaPsUGamys7Pxer3Q6uQ3WkFJcUlUzQRSBUBRmiHrwH46OtwAtHXomDTUzWARYs8efwMC2Ur+5L979+41LFOoqQKgKE1UU1PD4bwjdIrzf/I3a9A+VpKVlWVwMqUhdu3ahWbTIPbkN1oBGuzcudPIWCGlCoCiNNGhQ4fQdZ3O8b5T3+sS62Z/lI0jR6qMbRn4WvtAnPyGGUiCbdu3GZgqtFQBUJQm2r9/PwBdflQAOsf7KC45QWlpqUGplIY4fvw4+cfykW1+2sFVb6OzZ88eampqDEoWWqoAKEoT7d+/H5tZkGrXT32vthjs27fPqFhKA2RkZAAgU35aAGSKxOvxRk1jOFUAFKWJ9mVm0iXOgyZ++F5XVQAiwubNmxFWAUk/+0EKIPw/jwaqAChKE/h8PrKyDtAl3vuT78daJCmOH4aHlPAjpWTDxg3obfQfxv9rWYBk2LBxgxHRQk4VAEVpgiNHjuCqcZ/6xP9jXWLd7M+MnqmEkSYnJ4fiouJT/X9+Tk/T2Ze5j/Ly8hAnCz1VABSlCc50AbhWl3gfR47l43Q6Qx1LaYANG/yf7mXamQuATJNIKaNiGEgVAEVpgoMHD6IJTrV/+LHOcT6klOqGsDC1YcMGRIL4Yf7/zyWDiBGnCkVLpgqAojRBdnY27WIlljO8gzqcvDFMFYDwU1NTw9aMrfhSTy/cp2igp+isW78OKc98ltBSqAKgKE2Qcyib9nbPGX+WatexaP6xZiW8bN++HY/bU+f4fy3ZVlJcVNzin0NVAMJAdXU1jz76KPeNG8dnn31mdBzlLLxeL8eO5dP2DMM/AJqAtFjJ4cNqbYBws3HjRoQm/NM961FbINavXx+CVMZRBSAMHDx4kA0bNpC5dy9ffvml0XGUsygqKsLr8/3kBrCfS7F6yD8WXatLRYL1G9b77/41n+WBDv/6ABs3bgxJLqOoAhAGatsGSLONkhMnjA2jnFXtur9tbHUXgDY2nfx8tT5wOCkpKSH7YDZ6Wt3P24/5Un1kbMvA7XYHOZlxVAEIA7X9x3VHa4qLiw1Oo5xN7XOUZK37QJJklVRVO6Omp0wkqJ3WKVMbdmFXpkncNe4W3RZCFYAwUHtA8cUmU1Zais9XzwwFxXC1Z2wJlroPJAkx+k8eqxhv8+bNaFbthwVgzuZkW4hNmzYFM5ahVAEIA4WFhQirA2mNR9d1TqhhoLBWWVkJgKOeAhBrlj95rGKs2vYPvja+09s/1OVkWwhVAJSgKigoQLc4kDFxp75WwpfL5UITnPEegFpWkzz1WMV4eXl5FBUW1Xn3b130VJ29e/dSUVERpGTGUgUgDBw9dgxfTDy61V8A8vPzDU6k1Mfj8WDR6v8YadZ+eKxivNrZPGcqACJDIDLO/HzWtoXYunVrUPMZRRUAg/l8Pgry89GtcUhrPABHjx41OJVyVg0dRlDCwoYNGxBxAuJO/5koFYjSOp7Q1iAsLXc6qGEFQAhhE0JsEEJsE0LsEkI8a1QWIxUWFuL1epG2RDBZENZY8vLyjI6l1EPTNHS9/qGE2h9rmvqMZTS3283mLZvxpTVhcsXJthDfr/2+RbaFMPLVWQNcIaUcCAwCRgohLjQwjyFyc3MB0G2JAHitCeTk5BoZSTkLq9WKR//hIH8mbl2ceqxirIyMDGpcNch2TTuAy3aSwuOFLbK3k2EFQPrVTpGwnPzT8krsWRw6dAgA3Z7k/68tkexDh1rkp42Wwm63A+Cq5wOlyyt+8ljFON9++y3CLCC1af++tnCsXr06gKnCg6Hnp0IIkxAiAzgOrJRSntZ4QwgxTgixSQixqbCwMOQZg+3gwYOIGDtY/AcK3Z6Ms7qKlvj/2lLEx/uv1VR7674QUPuzuLgzDDorIeP1evl61df42vnA1MSN2IE2sPKLlS3ug5mhBUBK6ZNSDgI6AkOEEP3P8Jg3pZTpUsr0lJSzdHCKQAcOZOG1/XBniu5IBiArK8uoSMpZJCQkAFDprvvtU+ERP3msYoz169dTUV6B7Ny8A7feWSc3J5fMzMwAJQsPYXGFSkpZCqwCRhqbJLS8Xi8Hsw+iO1qf+l5tAVBryoavpKQkAMo9dZ8BlLsF8bEOzOazdR1Tguk/H/8HYRfQtnnbkZ0kwiz45JNPAhMsTBg5CyhFCJF08u924FdAVC2kmpOTg9fjwRf7QwHAHAP2RPbt22dcMKVeycn+Il1WU/fbp8ytnXqcYozc3Fw2rN+Ar5uv+Ue6GPB18rF8xfIWdae+kWcA7YCvhRDbgY34rwF8amCekNuzZw8AetxPh7a8jjbs2rXbiEhKA7Ru7S/Ype66zwBK3SaS27S8IctI8u6774IJZM/AjNvL3hKP28MHH3wQkO2FAyNnAW2XUg6WUg6QUvaXUv7JqCxG2b17N8JiQ1p/Ok7si0uhuLhIXQgOU3a7HbvNSml9ZwAe06lCoYTegQMHWLlyJb7uPgjUTNwE/7WAD/71AcePHw/QRo0VFtcAotW27TvwONqA+OknST3OP19t586dRsRSGiA5uRVl9VwELnMJNQRkEF3XWbBwAcSAPCews3Zkf4nX5+XV114N6HaNogqAQcrKyjicm4Men3baz3RHa4TJzI4dOwxIpjREcnIbyusYAnL5oMYnadWqoX2HlUD673//y/Zt2/H190FMgDceC76+Pr7+6mu+++67AG889FQBMMj27dsB8MWfYXqCZsIbm8LWjIzQhlIaLD4hgUrfmWf4VJ6cHZSYmBjKSAr+GysXv7IY0kB2C86cfdlXIpIEs+fMpqioKCj7CBVVAAyybds2hGY+7QJwLV98Ow5mZbXYNrSRLj4+vs4bwZzqJjBDVFVV8eSMJ/EID74LGtH3v7E08A71UlldydNPPx3RHV9VATDIps2b8calgHbmT5G+hHZIKdm2bVuIkykNYbPZ8PjOfISpOfl9m80WykhRzev18syzz5Cbm4t3qNd/924wJYDvfB87duzgpZdeitg7hFUBMEBpaSkHs7LwJbSv8zF6XCrCZD61jqkSXsxmM946lgT2yR8eowSfruvMnTuX9evWow/Wm9zzp7FkZ4neT+ezzz7jrbfeCs1OA0y9Qg2wZcsWAHwJHep+kGbCG5fGxo0tdzm6SCalPOsQg67XvWi8Ehi6rrNw4UI+//xz9HN1ZI/QfhKX/SS6S+e9997DarVy1113hXT/zaXOAAywceNGhNmKHtem3sd5EzuQm5vTYuYctyQulwtrHc3FrCffVTU1NaELFIV8Ph8vv/wyH330EXpvPeBTPhtEgDxPonfWeeutt3jnnXciajhIFYAQk1Kyfv0GPPFtQdT/69dPniG05EWpI1VFRQUO85k/4ceeXCy+vLw8lJGiitvt5plnnuG///0vel8dOeDsZ2RBI0AOkehddZYsWcKCBQvw+Zqw+IwBVAEIsZycHIqKCvEldjzrY3VHMiLGwYYNG0KQTGmMgvxjJMd4z/izJKuOAHXmFiTl5eU8MuURvvnmG/QBOvIXBh78awmQ6RK9j85HH33EU08/hcvlMjjU2akCEGK1B3Nf0tkLAELgTujA+g0bIuYTRTTQdZ28vDzS7Gd+TiwatLbD4cOHQ5ys5cvNzWXc/ePYsWMH+lAd2SeMhlsEyAESfZDO6m9X8+CEB8O+nYsqACG2dt06cLQ6tQD82fiSOlJVWcnu3ao5XLg4cuQI1U4XnePrLsqdY91k7t0TwlQt39q1a7lv3H0cKz6G9zJvs3v8B4vsJfFd7ONA9gHG3js2rFu6qAIQQk6nk4yMDDz1zf75GV9iBxCC9etPWyxNMcjWrVsB6J145iEggN5JXg7nHYn4O0XDga7rvPvuuzzxxBO4rC68V3ih/vkTxmsP3su9lHnLmDhxIh9//HFYXhxWBSCEtmzZgs/rxZvUqeH/yGxDj0tl7dp1wQumNMq6dWtJtkP72LqnefZP9p58rHremqOiooKp06by1ltv4evkwzvcC7FGp2qgRPBe4cWX4p+tNHv27LCbGaYKQAitX78eYbKgn6n/Tz28iR3Zv38fJSUlQUqmNFR5eTnr161jSIrr501cf6JLvI80h+TLL74IXbgWZt++fYy9dyxr161FH6wjh8jIu3MpBnzDfOj9dD7//HPuH38/eXl5Rqc6RRWAEJFSsub77/HEtwOtcatT+06eMahhION99tlneLw+hrV11/s4IWBYWxdbtm5RF4MbSUrJp59+yvg/jKegrADfcJ9/URejZ/o0lQB5rsQ3zEd2XjZj7x3Lt99+a3QqQBWAkMnJyaHw+PFTB/PG0B2tEdZYNZxgMI/Hw78++Cd9W3npmnD2WVlXdKjBLOD9998PQbqWweVyMXv2bObOnYs32Yv3V15oKevqtAPvlV6cNidPPvkkr732Gl5v3deRQkEVgBCpPXg3pQD8eDqo0S+YaPbJJ59wvLCI67s6G/T4RKvksvYu/ve/ZeosoAGOHDnC/ePv97d16KfjuySAq3mFi1jwDfeh99B5//33eeihhyguLjYsjioAIbLu1PTPprUI9iV2orqqSk0HNUhpaSl/fedt+iV7+UVyw4vwmG4uzELy2muvhuUskHCxdu1axo4dy6Ejh/AN8yHPjeAhn7MxnWwfMURnx+4d3DP2HsOmiqoCEAJOp5Nt27fjbsDdv3XxJbZX00EN9Oqrr1JVWcmdvavqvfj7c4lWyZiuVaxZ8z3ffPNN8AJGKCklf/vb33jiiSdwWp14r/RCO6NThYbsIvFe7qXUXcrEiRP59NNPQ55BFYAQ2Lp1Kz6vt0HtH+pktqLHpbJ+vWoLEWqrVq1i+fLlXNvVSae4xnf4HNW5hq4JOi+/9KK6L+BH3G43zz33HG+++Sa+jhE2xTNQkvzXBXxtfMydO5dXXnklpHf9qwIQAk2d/vlz3sSO7NuXSWlpaWCCKWd19OhR5r4whx6JOmO6Na23i0mDB86twFVVyZ/+9Ky6joN/Ou1DDz/EypUr/W2ch0bgFM9AqZ0q2lPngw8+YPqT00PWR8iwAiCE6CSE+FoIsUcIsUsIMdmoLMG2fv0GvPFpjZ7++XO+RNUdNJSqq6uZ+sTjSLeTB8+twNyMd0v7WJ3f9akkI2Mbr776auBCRqDjx4/zwIMPsGvXLvQLdWS/Fjze31AayMESfbDO92u+5+FHHg7JcrBGngF4gSlSynOAC4EHhRD9DMwTFPn5+Rw9egRvI9o/1EWPbYOw2FQBCAGv18uzzzxDTk4OE/uXk+po/uIul7R3M6qziw8//JB///vfAUgZefLz83ngwQc4fPQw3ku8yE7qwviPyZ4S3y997N6zm4mTJlJWVhbU/RlWAKSUx6SUW07+vQLYAzT/KBlmapd0bNb4fy2h4Ylry/oNG9WMkiCqXWJw7bp13N2niv6tAzdkc1svJ+eleFi4cAFfffVVwLYbCQoKCpgwcQKFJwrxXuoN2dKNEacjeC/ykn0om8kPTQ7quhJhcQ1ACNEVGAycNsVFCDFOCLFJCLEp3FurnsnmzZsRVgfSnhSQ7fkS21NcVBhWt5O3JFJK5s+fz+eff84N3Z1c2bH+O34bSxMwoX8lvZN8PDdzJqtXrw7o9sNVbQ//whOFeC/xQrLRicJc2x+KwBNTnwhaDyHDC4AQIg74EHhISnlaqZNSvimlTJdSpqekpIQ+YDNIKdm4aROeuHY0au5gPWrXEa5dV1gJHF3XmTdvHh9//DHXdnHxmyZe9D2bGBNMGVhOlzgPTz/1VIsvAl6vl+lPTicvLw/vL9XBv8Hagu8CHzt37GTW7FlBOes3tAAIISz4D/5/l1K2uEHRnJwcykpL8SUEbmKztCUgrLGnWhIrgeH1epk9e/apg/8tPZ2Bqtln5DDD44PL6BLnYcaMGaxYsSJ4OzPYm2++ybaMbfjSfWE37CMyBJQCpaCt0vxfhxHZSaL/Qufrr77mgw8+CPj2jZwFJIC3gT1SynlG5Qimbdu2AQS0ACAEnri2bNm6VV0HCBCXy8VTM2awfPlybuzuDPrBv5bDDE8MLqNvon8+/NKlS4O/0xDbvHkz77//PnoPHdkl/F6volQgPCf/FApEaXgVAADZRyLbS9544w2ysrICum0jzwAuBu4ErhBCZJz8c42BeQJu27ZtCGss0poQ0O364ttSeuKEug4QAKWlpTw0eTJr1qzhrj7V/KZ7/W2eA81uhkcHVXB+iptFixbx2muvoevNn3EUDmpqapjzwhxEvPAv2q40jQA9XUe36MyeMzugrw8jZwF9J6UUUsoBUspBJ//8z6g8wbB1awae2NSAjf/X8p28oWz79u0B3W60ycnJYfz949i/by+TBlQyopMxi3XEmGDygCqu6uji/fff5+kIWVD8bJYuXUpBfgHe87zRe5NXoFjBN8DHvsx9rFy5MmCbNfwicEtVUFBAcXHRqYN1IEl7EsJiC+u1RsPd5s2b+cP4+6ksyWfa4DIuSPUYmkcTcFcfJ7f3qubbb75l4oQHI7pthNPp5O//93dkOxl24/6RSnaW0Ar+uuSvAWsXoQpAkNQenPX4tMBvXAg8sSls374j8Ntu4aSULF26lClTppAkqng2vYxeSaHrvVIfIeCaLjU8PLCSQwcPcN+9YyO2++vKlSuprKhE79syhrPCggC9j87RI0cDtjaIKgBBsmfPHoRmRrcHZ86bHpfK4cO5IbldvKVwu93MnTuXRYsWMSi5hqfTS0mxh98B6rwUD8+kl6K5TjBx4gSWL19udKRG+3z554hE0XIWcwkTsoNE2ETAZo2pAhAku3btwhfbGrTg/Ip9cf57Ivbu3RuU7bc0RUVFTJo4kWXLlvHrrk4eGliJPUDj0u9l2nkv0x6YjZ3UKU7nT+ml9Ixz8fzzz7No0aKIaSJXWlrKzh078XXwqR4/gaaBr72P79d+H5DXgyoAQeD1etm3fz++2ODduKaf3HZmZmbQ9tFS7Nixg3vH3kPW/r1M+kUlN/V0oQXwwJRTYSKnonmN/s4kPkby+OAKru7kYunSpTzyyMOcOHEi4PsJtNrpzzJNzfwJBpkmqXHVBOS9rwpAEOTm5uJxu/HFtgneTsxWsCeqAlAPKSX/+c9/mDxpEhb3CZ5JL2VImrEXexvLrMGdfZyMP7eKXTu2c9+9Y8P+rC8rK8v/yb+V0UlaqJO/10DcE6AKQBDs27cP8C/m3lAxOWuJyVnbqP147cnszdzXqH8TLWpqapgzZw7z5s3j3FYu/pRe1qTFXMLFsHZunjq/DL2yiAcffIDPPvvM6Eh1ysvLQ4vVIPAnRQqAA4RJBGSdaVUAguDAgQMIzYy0Jzb432hVxWhVjVscWnckU5B/jKqqqsZGbNEKCwuZNHECn332Gb/p5mTKwEpiLZE/HNEtwcfMC0rpFe9i9uzZzJ8/PyyvC5SUlKBbI7fYhj0Bwi4CMhyoCkAQZGVloTuSQAT311t7hnHw4MGg7ieS7Ny5k/vuHUv2gX1MHlDJjT0CO95vtPgYyeODKhjV2cVHH33ElEceCbsV4qqd1Uhz5BfccCbNEqfT2eztqAIQBAeyDuK1BX8AVHf495GdnR30fUWCzz//nMmTJmKpKeWZ9FLDb+4KFpMGd/R2cn+/KnZuz2DcffeG1YcAqasVvkIhEL3AVAEIsPLycspKT5w6OAeTjIlDmCxRXwB0XefPf/4zs2bNoldCDc9eUErHCB7vb6hL2rt58vxyXKXHeeAP4wN2c1Bz2ew2/3p/StAIn8BmszV7O6oABFhOTg4A0pYU/J0JgW5PPLXPaFRTU8MzzzzD3//+dy7vUMMfB1UQ1wLG+xuqR6KPZy8oJcVczROPP87HH39sdCQS4hPQPBFyaPGA3W7nxhtvxG63Q6ScNNZAfHx8szcTIc9S5Kg9GOsBWgHsbHzWRLIPRWcBqKioYMojj7Bq1Spu61XNPX2rm7Vwe6RqbZPMOL+MAa3dvPzyy/zlL38xtFV4SkoKOIFIqMMeGD16NJMmTWL06NGRUQC8oLt1ArFAlurRF2CHDx8GzYS0xoVkf7o9ieK8LFwuV0BOCSNFSUkJUx55hJxDB5nQv5IL20bCOzd4bGZ4eEAlSzIdvPfee5SVlfHwww9jMoV+LmaHDh2Qbgk1QLi/JC2wbNky4OR/rQbnaYhK/386dmz+OuNR+HkpuPLy8sCWEPQZQLV0m3+tgSNHjoRkf+GgsLCQiQ8+QF5ONlMGVkT9wb+WSYN7+lZzXVcnn3zyCbNmPR+wrpGN0a1bN/9fykK+68az+DuXLl261D+rxmJ0oLOrXbSma9euzd6WKgABlnv4MN6Y5o/NNZS0+e81iJbFYYqLi3lo8iSKCo7x+OAyftFaXW38MSHglp4ubu7hZOXKL5gzZ07IF5jp3bu3P8sJNRUoKE5AjDWGTp06NXtTqgAEkJSSo0ePIm2BXQGsPro1es4AKisrmfLIwxzPP8qjg8roHQZtnN/LtJ/qBfTcpriAN4Vrquu7ubixu5Ply5ezcOHCkF4TSEhIoGOnjohCVQCCwVRsov+5/QMyvKcKQAAVFxfjcbtPDcuEhDkGEWNv8QXA4/Ewffo0cnIO8fCAcvqEwcEf/I3gnD4Np09jb6klKE3hmuo33V1c08V/w9g//vGPkO47/fx0tCINwuNpajlcIE9IBg8eHJDNqQIQQLXDMIFeA/hsfDHx5LXwArB48WK2bs3gvnOq6J+shn0a6taeTi5Mc/PGG2+E9D6BoUOHIr0SIndRs7Ak8v1nVRdeeGFAtqcKQADVFoCQngHgHwbKzW1+Y6hw9cUXX/Cf//yHazq7GNbObXSciKIJuK9fFZ3jdWb+6VkKCgpCst/09HSsNisiTw0DBZLIE7RJaXPqOktzqQIQQHl5eSC0kE0BraXbEiguKqSmxphFzYOpqKiIeS+/RK8kHzf3bH7vk2hkNcGkX1TgcVUz94UXQnI9wGq1MuziYZiOmNQwUKDUgCgQ/OrKXyFEYAqrKgABdPjwYbCHbgpoLd3ecmcCLVq0iBqXk3HnVEblTV6B0tahc0uPKjZu2hSw5QTPZtSoUcgaiTiqzgICQeQI0GHkyJEB26Z6SwXQoZxcvNaGt4AOlNqpoLm5uSHfdzDt3LmTVatWcV2XatrFtvzePsF2Zccauifq/OXNP4fkbPH888+nbbu2aAfUYabZJJgOmji3/7l07949YJs19JkRQrwjhDguhNhpZI5A8Hg8HDmSF7IWED+mn+w7dOjQoZDvO5iW/PWvJFrhms4uo6O0CJqAW3tUcbyw6NTdr8FkMpm48YYb/ReCG7fUhfJzR0FWSP/vM4CMLs1LgMCdzxgoNzcX3edDtxuwDp7JjLAnhlVL4ObKzc1lw8aNXNXBiU01LAmYfsleuifq/PvDpSG5FnDttdcSGxeLtsfoQ00Ek2Daa6Jtu7ZcdtllAd20oc+KlPJboMTIDIGyf/9+oHHLQAaSp4UtD7lixQo0AZd3aHkXto12ZQcnuYfzQrK2sMPh4NZbbkUcE+osoKmOAiXwu7t/h9kc2E9DdRYAIUSd9xkLIS4JaIp6CCHGCSE2CSE2FRYWhmq3jZaZmYkwWRq1DGQg6Y7WFOQfo7y83JD9B9p3q7+ld5KXRGsktJSMLOeneNAEfPfddyHZ30033URCUgKm7abI6BAaTnQw7zDTsWNHRowYEfDN13cG8I0Q4o9CiFMlRwiRJoT4GzAv4EnqIKV8U0qZLqVMD0T702DZvmMn3tg2IZ8BVMsXnwrArl27DNl/IFVUVHAw+xC/SFZN3oIhziLpluBj+7ZtIdmfw+Hg/vvuhyIQh8NrRpBMkkjLyT8pEpkUXhVKHBDICsnEiRMD/ukf6i8A5wM9gK1CiCuEEJOBDcBaYGjAk0SwyspKsg7sxxeXZlgGPTYVNI2MjAzDMgTKvn3+oazuCeF/x6/TK36yoIjTG14HuLp0j/eQmZkZsh5B11xzDb169/KfBYRRXZeDJCQBSaAP1/1fh4tqMO0yMfTCoQG78/fn6iwAUsoTUsr7gbeAL4DHgIullK9KKdWcvB/JyMhA13V8ie2NC2Eyo8emsmHjRuMyBMixY8cA/9z1cFftFT9ZUKQ6QgpAW4eOq6aGsrLQ9Gw2mUw89uhj4AKxLTJ+R4aSoG3WsJgsPPzQwwG78evn6rsGkCSE+DPwe/wzdZYCnwkhrgjUzoUQ/8B/RtFHCJEnhBgbqG2H0tq1axHmGHQDzwAAvIkdyTpwgHC+VtIQxcX+q4WtrOFfABxmybJly1i0aBHLli3DYQ6jT5D1qP3d1v6uQ6Fv377cdtttaNkaHAvZbiOSOCQQ+YI/jP8D7dsH74NlfUNAW4D9QLqUcoWU8iHgTuC5kwfuZpNS3ialbCeltEgpO0op3w7EdkPJ6/Xy7erVeBI6gGZsJ0hvqy4ArF692tAczVVTU4NJEBF3/trN8icLitgjpABYTf6cLldo77H4/e9/T5euXTBvMoO6vePMKsCUYWLQ4EGMGTMmqLuq7y12qZTyJSnlqYFYKWWGlPIi4KugpoogW7Zsoay0FG/rwN2d11TSngSOVqxYudLoKM2i6zqoUYKgqh1RCPViMVarlWefeRaTz4Rpg5oVdBofmNebcdgcPDn9STQtuJ+C6rsGUGdjGSnlX4ITJ/L873//Q1is+JKavzpPswmBu3VPdu/aFdFtIWJiYvDpoKuDQ9B4dH8FiImJCfm+u3fvzkOTH4ICEHtUpf8xkSGQJyQznpxBampq0PcXASfZ4au4uJhvvvkGd+teoIXH7aqelF6gaXz00UdGR2my+Hj/kpqRckE1ElV5/L/b2t91qF133XVcffXVaLs0/41OCuKgQDuocccdd3DRRReFZJ+qADTDv//9b3w+HU/aOUZH+YHFgbdVdz5dtixibwpr1crfTqOsRhWAYClz+3+3SUlJhuxfCMGjjz5Kz149MW80Q2S+VAOnEExbTaSnp3PvvfeGbLeqADRRRUUFS5d+iDe566lunOHC3X4ANS4XS5cuNTpKk9Se+ha51MszWIpcGvFxDhwOh2EZrFYrs2fNJsGRgHmNGaK160clmNeZad++Pc8++2xA1vptKPUOa6K///3vOJ1OPO0HGR3lNNKRjDe5K+//85+cOHHC6DiNVjvt7bgzfNbXbWmOV5to176D0TFIS0vjhTkvYK4xY1oThYvHuMG8xkysJZa5L8wN+ZCcKgBNcOzYMf71r6V42vRAjzWm+dvZuDum43K5eOedd4yO0mitW7fGYbdxtEq9PIPlmNNC585djI4BQL9+/ZgxYwaiRKCt16JnZpAPTN+b0Ko1Zs+aTadOoZ9Iot5hTbBw0SK8usTTMd3oKHWS9iQ8qf345JNPyMzMNDpOowgh6NqtG3lV4XFhvaWp9kKRk4AuLNJcw4cPZ8KECYgjArFVtPwiIPEXuyKY8eQMBg4caEgMVQAaadWqVXy/Zg017QeHfO3fxnJ3PB8sdua88AJeb/j31fmxHj16kltpIUStaqJKboW/sPbo0cPgJD910003+e8UztIQu1vwBAAJYrNAHBFMnDCRK64IWHOFRlMFoBFKS0t56aWXkbFt8LTtH7DtxuSsRasuRqsuxrb7U2Jy1gZmw+YYnF1+SdaBA/ztb38LzDZDpE+fPlR5JMed6iUaaNnl/msrvXv3NjjJ6caPH8+oUaPQdmuI/S2zCIidAi1b48477+Smm24yNIt6dzWQlJJZs2dTXlmJs/ulEMA79LSqYoTPg/B5MFXko1UFrj+LL7kb3tY9WLJkSUS1iu7Xrx8A+8vUMFCgHSg3k5rSmtatw+/6lRCCxx57jGHDhqFlaIhDLasIiL0Cba/G9ddfH9LpnnVRBaCB/vWvf7Fu7VpqOg1BOpKNjtMoNV0vQo+J5elnnqGiosLoOA3SrVs3HHY7maWqAASSlJBZFsOAgYONjlIns9nMM888w3nnn4e2SYM6exJEFnFAoO3QuPLKK3n44eB1+GwMVQAaYPv27bz2+ut4W3XBm9bP6DiNZ7bi7HE5xwsLmfnccyHv/9IUJpOJgYMGsrvUanSUFuVotUapCwYNGmR0lHrFxMQwe9Zs+vXrh3m9OeK7h4pDAm2rxkUXX8T06dNDOte/PqoAnEVhYSEzZjyFHhNHTfdLf+iiFWH0uFRqOl3IurVrWbJkidFxGuSCC4ZQUAUF1eplGig7ii0ApKeH7wy2Wna7nRfnvkiPHj0wrzXDcaMTNY04LNA2aZx3/nk8+8yzQVnZq6nUO6seLpeLJ6ZOo7SiguqeV4I5sj+NetPOwdOmF0uWLOHrr782Os5Z1fZD2VxoMThJy7GpMIZuXbsEtcd8IMXHxzPv5Xl07tTZf7dwkdGJGukoaBs0+vfvz+xZs7Faw+sYogpAHXRdZ9asWezftw9n9+ERN+5/RkLg7jYMPT6N555/nj179hidqF7t27enV88erCuwGR2lRSh2CfaVmrls+OVGR2mUpKQk5s+bT7u0dv4iECk3txeAea2Z3r168+LcF7Hb7UYnOo0qAHV48803WbVqFTWdL8DXKjzumAwIzYSz16/waDb++Pjjp5ZfDFdXjxzFwXKNw5Xqpdpc3x2zoku4+uqrjY7SaG3atGHhgoW0SWqDebUZQrOSZdMVgvl7M127dmXey/OIjY01OtEZqXfVGXz88cf83//9H57Uvnjb/sLoOIFnsVPd+yrKq5w8+thjYT0zaMSIEcRYLKw4rM4CmsOrw5dH7Zw3eDAdOhjfA6gp0tLSWLRwEUmxSf4iUGl0ojqU+Pv7tG/bngXzF5CQkGB0ojqpAvAza9asYd68efiSOuPuelHEXvQ9G2lvRXXPX5GXd4QnnphKTU14tmJMSkpi5KhRfHfMygnVHrrJ1ubHUOKEW2691egozdKhQwcWLlhInCUO87dmqDY60c+Ugfk7M6mtU1m0cNGp1ubhShWAH9m1axdPPf00emwbXD0vB9Gyfz16Qjuc3S9lx47tPBfG00Nvv/12pDDx0cHwG0ONBB4dPsyOpXevngwdOtToOM3WtWtX5s+bj03a/EUgXNYWrgTLagtJsUksXLCQlJQUoxOdVcs+wjXC4cOHeeyPj+Mx2ajuPQJM0THzxNe6BzWdh/DNN9+wePFiZBg232nfvj2//s1vWHXUyqHy8Jg/HUn+l2OjyAn3j/9D0NeYDZXevf0XVs01ZszfmcFjcCAnmFebcZgdLJi/IGJmWbWMV0MznThxgkemPEqVy0N176vBEl2fNL1tf4Gn7bl8+OGHfPDBB0bHOaN77rmHxMRE3t4bhzc8T1TC0rEqjf9k27nsssu44IILjI4TUAMGDOC5mc8hygWm7w1cS8DtH/aJ8cbw0osv0a1bN4OCNF7UFwCXy8UfH3+c44WFVPe+KuxW9woJIXB3vhBvq668+uqrrFq1yuhEp4mPj+fhR6aQXa7x4cHwuSDcJd6H3aRjN+n0TfLQJT58VjTx6PDqrnhsjlgmT55sdJyg+OUvf8m0qdPgOIgNTWsjLZMkMqmJZ74+MK01oVVoPP/c86d6WEUKQwuAEGKkECJTCHFACPFEqPev6zrPPfc8mXszcXYfjh6XGuoI4UMIanoOR49PY+bM59i9e7fRiU4zfPhwRo8ezaeH7GQUhcfdlHf2cdIl3keXeB9PpldyZx+n0ZFO+ds+O4fKNaZOm06bNm2MjhM0I0aM4IEHHkDL0xDbGz9RQA6SyEFNqRwgNgo4DtOmTWPIkCGN34bBDCsAQggT8CowCugH3CaECGn5fOedd/j222+o6TwEX3LXUO46PGlm/z0CJhuPPzGV48fD7977yZMn06NHd17dmaDuDajHilwrX+bZuOWWWxg2bJjRcYLulltu4YYbbkDbpyEOhGa2mNgp0A5rjB8/nquuuiok+ww0I99BQ4ADUsqDUko38D7w61DtfNWqVbz77rt4UnrjDWBv/4hnsVPd6yrKKyuZNm162E0PtdlszHlhLo6EJF7MSKRQrRdwmg0FFv6238HFF1/E+PHjjY4TEkIIJkyYwEUXX4SWoUF+kPd3yN/W+brrruO2224L7s6CyMh3Twfg8I++zjv5vZ8QQowTQmwSQmwqLCwMyI5zc3N5ftYs9PhU3F0vbrFz/ZtKOlrh7DacffsyWbx4sdFxTpOamspLL8/DbXIwe2siJS71/NXaWmTm1Z1x9Dv3XJ566umw6ToZCiaTiadmPEX37t39HUSDdX9jMZi2mDjvvPPCpq1zUxlZAM70WzttIE5K+aaUMl1KmR6IebU1NTXMmPEUbh+4elwJWvS8QRrDl9wFd7sBfPLJJ3z55ZdGxzlNjx49eOnleVRKGzO3JKmVw/B/8l+wPZ4evXoxN0x7zwSbw+Fgzuw5xNniMH8fhOmhTrCstZCWmsaf/vSnsOrs2RRGvmvygE4/+rojcDTYO3377bfJzj6Is/tlSGt49ucIF56O6ejxabz40ksE6uwrkPr168e8+QtwabE8tzkxqq8JfHM0hld2xtG3bz/mz19AXFx4r1cdTG3btvVPD630t2EO2ALzOpjWmzDrZmbPmh3WLR4aysh3zEaglxCimxAiBrgV+CSYO9y1axfv//OfeFL74kvqdPZ/EO00DVf3S3G63Myd+2JY3iTWr18/Fi1+BeFoxczNSewuiexPZI0lJXyYZeMvu2M57/zzeXnePOLj442OZbjBgwczbtw4RJ5AZAVmiEbsElAIf3zsj/To0SMg2zSaYQVASukFJgDLgT3AB1LKoC1a6/P5eHnePESMA3fnMJuu5XNjt9u58cYb/aftPrfRiU6RtkRcHc5j/fp1rFmzxug4Z9SjRw/e+PObpLbvxAsZ8aw6EmN0pJBw++DPuxx8lG1n5MiRzJ37Ig6Hw+hYYePWW29lyNAhmLabmt899DhoezVGjx7NiBEjApIvHBh6ziyl/J+UsreUsoeU8vlg7mvFihUc2L8fZ6chYAqvA4Twuhk9ejSTJk1i9OjRCG/4FAAAb9q54GjFosWv4PV6jY5zRmlpabz62uucd146b+2J5W+Zdnwt+I7hEzWC57ck8F2+lbFjxzJ16tSIH48ONE3TmD5tOvHx8Zg3mqGprwcPmDeZad+hPZMmTQpoRqNFxaCp1+vlr0uWIONS8CV3NzrOaaQ5hmXLlrFo0SKWLVuGNIdXgULTcHVMJ//YUVauXGl0mjrFx8fzwty53HjjjXx+2MYLGfGUuyN3hkZd9peamLExiSMuOzNnzuTuu++O6JkowdSqVSsef+xx5AmJ2Nu035HYLqAaZjw5o8VdWI+KArBmzRryjx2jpt3A8JzyaYrB6XSydOlSnE5n2J2hAPiSOiNjW/OP998Py2sBtcxmM5MmTWLq1Knsr7AxY2MSWWUtY6aXlLDysJXntiTgaNWW19/4M5dddpnRscLeJZdcwhVXXIFprwnKG/mPi0A7qHHzzTdz7rnnBiWfkaKiAHy6bBnCGouvVWejo0QuIXCn9OVQdjaZmZlGpzmrUaNG8dprr2OOT2Hm5gS+yIshjOvWWbl88PouB/9/poMLhgzlL2+93WIuRIbCpEmTsFltmDJMDZ8VpIN5q5k2KW245557gprPKC2+AFRXV7Nx40Zqkru3+P7+weZt3R2EYPXq1UZHaZA+ffrw1tvvcP4FQ1iyN5bXdzlwhU+vtgY7WqXx9MYk1hZYuffee5kz5wU106eRkpOTuXfsvVAANHAVVHFIIEslEydMbHFDP7Va/BFx+/bt6D4fvsSORkeJfGYrelwqGzduMjpJgyUmJvLCC3MZO3YsawusPL0xiSNVkfOyX5dvYcbGJCq1BF5+eR533XVXi+npH2pjxoyhXft2mHeaz34W4AXTbhPn9j+X4cOHhyKeIVr8KykrKwsAPTb8V+eJBL7YFLIOZoX1dYCf0zSNu+++m5dfnkeVKZGnNiaxriBwC/7UdgMNJK8O72baeWVnHL36nMPb7/yV9PT0gO4j2pjNZsbdNw5ZJv23odZDHBRIp+T+cfe36AvsLb4AHD16FBFjh3CbWROhdGsCHrebkpISo6M0Wnp6Om+/81d69j6HV3bE8V6mPSCLy9zZxxnQNtAlLv8UzxWHbdx8880sWvwKqalR3Ko8gIYPH07Hjh0xZ9ZzFqCDab+JAQMHMGjQoFDGC7kWXwBqamqiZnnHkDD555qHW5fQhkpJSWHR4sXccMMNLD9sY87WBMrCaLH5vSfMzNiYRJ7LwbPPPsuECRPU/P4AMplM/jWmT0ioo7uJyBPIaskdt98R2nAGaPEFQAhBRE//UALOYrEwefJknnzySbKr/FNFD4bBWsMrD1uZvSWehDbt+fObf+Hyyy83OlKLdNVVVxEbF4s4eObCrx3UaNe+HUOHDg1xstBr8QWgdevWyJoqVQQCRLirAP/vNdKNGDGC115/A0tCCjM3J7I235gzRa8Of93jn+I55MIL+fNf3qJr166GZIkGVquVkVePxHTEBD+/6b4SKITrrr0uKi62t/j/w44dO4LUEa7mNgNRADTnCVq3boPVajU6SkD06tWLN//yFuf068+rO+P46KAtpJ8Vqr3wUkY8Xx6xcvvttzNr1uyo7uQZKiNGjEDqEnHkp2cBIlec+nk0aPEFYMCAAQCYKoK8RFA0kBJLZQGDBw8yOklAtWrVivkLFnD11Vfz4UE7b+1xhKSPULFL8KdNiewpi2Hq1KmMHz8+qhZwMVLfvn1Ja5t2WgEwHfVP/YyWi+4tvgB06tSJtLS2mEuyjY4S8bTKAmRNVUQufn02FouFadOmcffdd/PNUSuLdsThDuJNY/nVGjM3J1Hic/Diiy8xatSo4O1MOY0QgmEXD0M7rkHt8+wEeUIy7OKWv4ZyrRZfAIQQjBo1ElPZEYQrWGvERQfz8UysVhuXXnqp0VGCQgjB2LFjmTx5MpsLLby8LZ6aIBSBvEqNmZsT8VjiWbhosZrfb5AhQ4YgfRKK/F+L4+LU96NFiy8AANdeey0msxnLsW1GR4lYoqYCS3EW11wzqsX3nL/hhhuYNm0au0+Ymb8tPqBnAkeqNGZtTcQc24pXXn2NPn36BG7jSqMMHDgQIQSi8OQwUCHExsVGVY+lqCgAqampXHfttViK9iGcpUbHiUiWvC2YTBp33NHy50YDjBw5kieemMquE2YW74gLyDWBIqfG7K2JmOyJLFy0mC5dujR/o0qTORwOunbriijxFwBTiYlz+50bFbN/akXN/+nvfvc77DY71py1akpoI2nl+ViK9nP7bbdFzcUx8HcUffjhR9haZGFJpqNZL5tKj2DutgS8JgfzFyykc2fVmTYcnNP3HExlJvCBLJf07dvX6EghFTUFIDk5mfH3j8NUdgRzYXi1M9ZjWyNNFqTJgi++LXpsGM2x93mwH1pNaload955p9FpQu43v/kNd955J18fsbL8cNOmvvp0eGVHHIUuM7Nmz6F79/BblCha9ejRA92l++8KlkTV8A9EUQEA+PWvf83556djy12HqD5hdJxT3F1+ie5oje5ojavftbi7/NLoSH5SEnPoe4SrnBlPPonNZjM6kSHGjh3LsGHD+L/9DvacaHxbhqUHbewsMfPIlEdbfG+ZSFM7DGdeb/7J19EiqgqApmlMnz6N+Lg4HAe+AG9k9rMJFfPxPViK9nPXXXcxcOBAo+MYxv+6mU6H9h14fXc8VZ6G9w7aWWLmv4fsXHvttYwePTqIKZWmGDhwIDfeeCOjfjWKO+64I+ruwBaR1NY3PT1dbtrU/F7027ZtY/JDD+GJa4ur99UQBhd9bLs/BcDV71qDk/hpZUewZy7nwqFDmDNnTlRdGKvL3r17+cMfxnNhipM/9K8+6+OrvfDE+lbEtu7A2+/8NWrPoBTjCSE2SylPm28cle/qgQMH8uiUKZjKjhCTvVpdFP4ZraoYx4Ev6dKlM0899ZQ6+J/Ut29f7rzzLtbkW9lVcvahoH9n2TnhgulPzlAHfyUsGfLOFkLcJITYJYTQhRCG3AVz7bXX8vvf/x5L0X5icjeoInCScJXh2Lec5KQEXn7pJdWX5mfuuOMO2rdry7v74tDreckcqdJYkWfjuuuup1+/fqELqCiNYNRHu53Ab4FvDdo/4J8aOmbMGCz5O7DkbTYySlgQrgpi935GvM3C/HnzomrKZ0NZrVbuH/8HjlQK1hyre5GhD7PsWK3+NXwVJVwZUgCklHuklIbPxRRCMHnyZEaPHk3M0QwshzdF7ZmAcJXhyPwfDgssWDA/6i6GNcZll11Gzx7d+STHccazgCNVGhuOx3DTzbeQlJQU8nyK0lBhP7grhBgnhNgkhNhUWFjHEj7NoGkajz322I+KQPQNBwlnKbF7/0ecRbBo4UJ69uxpdKSwpmkat91+B8eqBNuKT78WsDzXhsVi5oYbbjAgnaI0XNAKgBDiCyHEzjP8+XVjtiOlfFNKmS6lTE9JCc7C7rVFYMyYMcQc20FM9ncgQ9APOAxoVUXE7llGvM3CK4sX0atXL6MjRYTLL7+cVkmJrDry05vDXD74vsDGlVf+ilatWhmUTlEaJmiLjUopfxWsbQeDpmk89NBDxMXF8d577yF8bmp6DActNP3Zjbj7Vys/imP/F7RJTmLB/Pl06tQp5Bkildls5uqRo/jXP9+n0lNNnMV/1ril0ILLK7nmmmsMTqgoZxf2Q0ChJITgvvvu48EHH8Rcko0tczl4f75mXHC4u/wypHcAm4qzcWQup2P7trz+2mvq4N8El19+OT4JGUU/LCW56XgMrZNbnVqISFHCmVHTQMcIIfKAXwLLhBDLjchRl1tuuYXp06djqSrAsXcZwn32m34iiTl/N7YDX3HOOX15/bXX1GyfJurTpw+tk1ux9WQB8Oqwo8TKRRcPU/dOKBHBqFlAH0kpO0oprVLKNCnl1UbkqM/VV1/NC3PmYPNW4djz35bRRlpKLLkbseZ8z0UXX8SC+fNJSEgwOlXE0jSN89MvYE9pDFLCoQoTTq9UC7woEUN9TKnH0KFDeeWVxSRYTcTu+RQtktcV1n1Ys1YRc2wb119/Pc/NnKnuTg2AgQMHUl4DBU6NfaX+S2pq+EeJFKoAnEWfPn14889v0C61NY69n2GKxLWFvW5s+5ZjLs7i3nvvZcqUKZjNQbv+H1XOOeccAA6WmzhYbia1TWtatw6jdt6KUg9VABqgffv2vPH66/Tt2wfb/q8w5+8yOlKDCXcVjr3LsFQWMG3aNO666y6EaHg3S6V+Xbp0wWTS2F5kIbvCQs/eaolHJXKoAtBASUlJLFywgIsuvghrzloshzeG/Q1jwlmKY8+n2HxVzH3hBUaOHGl0pBbHYrHQvVs3vsu3UlAt6N27t9GRFKXBorIddHN4vV4WLFjAJ598gielN+5uw0CEXx3VKo/j2LeCeIeNl196US0+HkQlJSUcPnwYTdPo3bs3VmvTVg5TlGCpqx20GghuJLPZzJQpU0hOTmbJkiVoXheunleAFj6/SlNpHvYDX5Ka0oYF8+fRoUMHoyO1aMnJySQnJxsdQ1EaLfw+ukYAIQT33HMPDz/8MKbSw9gzV4AvNDeMnY2pJBvb/hV069KZN15/TR38FUWpkyoAzTBmzBienD4dc2U+jr2fG77EpKloP7YDX9HvnHNYvHiRmo2iKEq9VAFophEjRjBz5kzMzmIcez8Dj8uQHObjmdiyvmXwoEHMnzeP+Ph4Q3IoihI5VAEIgEsuuYTZs2djcZfhyAz9mYD5eCbW7NVccEE6c+fOxW63h3T/iqJEJlUAAuTCCy9k9qxZmGtKTxaB0FwTMBUdwJr9HenpFzBr1iw1A0VRlAZTBSCAhg4dyvPPPYfJWYJ9/0rQvUHdn+lEDraD3zBo0CBmz1YHf0VRGkcVgAC76KKLmDZ1Klr5MawHvg7awjJaRQH2A1/Tu1dv5syZrQ7+iqI0mioAQTBixAgmTJiA+UQOMbkbAr594SrHceAL2rZN5aWXXsThcAR8H4qitHyqAATJzTffzG9/+1ss+TsxH98buA173Tj2r8ARY+alF19Ui44ritJkqgAE0cSJEzk/PR1bzlq0yuPN36CUWA9+g+YqZ9bzz6lVvBRFaRZVAILIZDLxzNNPk5LSBvuBr5o9PdScvxPziRweeOABBg8eHKCUiqJEK1UAgiwxMZE/PfssmteJNfu7JncQ1aqKsOZtZNiwYdx0000BTqkoSjRSBSAE+vXrx71jx2IuycZUcrDxG9B92LO/pVVSKx5//HHVz19RlIBQBSBEbr31Vnr36YM9Z12j20VYjm6DqhIe/+NjJCYmBimhoijRRhWAEDGbzUx94gmEr4aYvM0N/neipgLrse1cfvkVXHTRRUFMqChKtFEFIIR69OjBmDFjsBTuRVSXNOjfxORuwGIx8eCDDwQ5naIo0UYVgBD7/e9/j91mb9BZgFZVhLkkm9tvu43U1NQQpFMUJZoYUgCEEC8KIfYKIbYLIT4SQiQZkcMICQkJ3H77bZhP5KBVFdX7WMuRrcTGxXHzzTeHKJ2iKNHEqDOAlUB/KeUAYB8w1aAchvjtb3+L1WbDcmxHnY8RzlLMJ3K46cYbiYuLC2E6RVGihSEFQEq5QkpZ2ypzHdDRiBxGiY+P59fXX4+5JBvhrj7jYywFuzGbLYwZMybE6RRFiRbhcA3gHuCzun4ohBgnhNgkhNhUWFgYwljBdf3114PUMRftP/2HupeYkiyGD7+MVq1ahT6coihRIWgFQAjxhRBi5xn+/PpHj5kOeIG/17UdKeWbUsp0KWV6SkpKsOKGXOfOnTm3f39iig+c9jPTiVykp4bRo0cbkExRlGhhDtaGpZS/qu/nQoi7gWuBK6VsYn+ECPerK69k186FiOoTSMcPn/TNJdkkJiUxaNAg48IpitLiGTULaCTwOHC9lPLMg+BR4NJLLwXAXJr7wzd1H5byPC695BJMJpNByRRFiQZGXQN4BYgHVgohMoQQbxiUw1ApKSl079ETc1neqe9pFQVIr4cLL7zQwGSKokSDoA0B1UdK2dOI/YajC9LPJ/tfSzEV+5vEmU/kIIRQ7Z4VRQm6cJgFFNWGDBmC1H3YDnyF7cBXmIuzOKdfPzX3X1GUoDPkDED5QXp6Ov/4xz9wu92nvpeWlmZgIkVRooUqAAYTQtChQwejYyiKEoXUEJCiKEqUUgVAURQlSqkCoCiKEqVUAVAURYlSqgAoiqJEKVUAFEVRopQqAIqiKFFKRFIjTiFEIZBjdI4gagPUv06kEq7UcxfZWvrz10VKeVo//YgqAC2dEGKTlDLd6BxK46nnLrJF6/OnhoAURVGilCoAiqIoUUoVgPDyptEBlCZTz11ki8rnT10DUBRFiVLqDEBRFCVKqQKgKIoSpVQBCANCiJFCiEwhxAEhxBNG51EaTgjxjhDiuBBip9FZlMYTQnQSQnwthNgjhNglhJhsdKZQUtcADCaEMAH7gKuAPGAjcJuUcrehwZQGEUJcClQC70op+xudR2kcIUQ7oJ2UcosQIh7YDPwmWt5/6gzAeEOAA1LKg1JKN/A+8GuDMykNJKX8FigxOofSNFLKY1LKLSf/XgHsAaJmiT5VAIzXATj8o6/ziKIXoKKECyFEV2AwsN7gKCGjCoDxxBm+p8blFCWEhBBxwIfAQ1LKcqPzhIoqAMbLAzr96OuOwFGDsihK1BFCWPAf/P8upfy30XlCSRUA420EegkhugkhYoBbgU8MzqQoUUEIIYC3gT1SynlG5wk1VQAMJqX0AhOA5fgvQH0gpdxlbCqloYQQ/wDWAn2EEHlCiLFGZ1Ia5WLgTuAKIUTGyT/XGB0qVNQ0UEVRlCilzgAURVGilCoAiqIoUUoVAEVRlCilCoCiKEqUUgVAURQlSqkCoChNJPy+E0KM+tH3bhZCfG5kLkVpKDUNVFGaQQjRH/gX/h4yJiADGCmlzDIyl6I0hCoAitJMQoi5QBUQC1RIKWcaHElRGkQVAEVpJiFELLAFcAPpUsoagyMpSoOYjQ6gKJFOSlklhPgnUKkO/kokUReBFSUw9JN/FCViqAKgKIoSpVQBUBRFiVLqIrCiKEqUUmcAiqIoUUoVAEVRlCilCoCiKEqUUgVAURQlSqkCoCiKEqVUAVAURYlSqgAoiqJEqf8Hx3Ih0akOp3UAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.violinplot(x=\"Y\", y=\"X\", data=df);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Improper Prior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will model the outcomes Y as coming from an OrderedLogistic distribution, conditional on X. The `OrderedLogistic` distribution in numpyro requires ordered cutpoints. We can use the `ImproperUnifrom` distribution to introduce a parameter with an arbitrary support that is otherwise completely uninformative, and then add an `ordered_vector` constraint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [00:03<00:00, 258.56it/s, 7 steps of size 5.02e-01. acc. prob=0.94]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "   b_X_eta      1.44      0.37      1.42      0.83      2.05    349.38      1.01\n",
      "    c_y[0]     -0.10      0.38     -0.10     -0.71      0.51    365.63      1.00\n",
      "    c_y[1]      2.15      0.49      2.13      1.38      2.99    376.45      1.01\n",
      "\n",
      "Number of divergences: 0\n"
     ]
    }
   ],
   "source": [
    "def model1(X, Y, nclasses=3):\n",
    "    b_X_eta = sample(\"b_X_eta\", Normal(0, 5))\n",
    "    c_y = sample(\n",
    "        \"c_y\",\n",
    "        ImproperUniform(\n",
    "            support=constraints.ordered_vector,\n",
    "            batch_shape=(),\n",
    "            event_shape=(nclasses - 1,),\n",
    "        ),\n",
    "    )\n",
    "    with numpyro.plate(\"obs\", X.shape[0]):\n",
    "        eta = X * b_X_eta\n",
    "        sample(\"Y\", OrderedLogistic(eta, c_y), obs=Y)\n",
    "\n",
    "\n",
    "mcmc_key = random.PRNGKey(1234)\n",
    "kernel = NUTS(model1)\n",
    "mcmc = MCMC(kernel, num_warmup=250, num_samples=750)\n",
    "mcmc.run(mcmc_key, X, Y, nclasses)\n",
    "mcmc.print_summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `ImproperUniform` distribution allows us to use parameters with constraints on their domain, without adding any additional information e.g. about the location or scale of the prior distribution on that parameter.\n",
    "\n",
    "If we want to incorporate such information, for instance that the values of the cut-points should not be too far from zero, we can add an additional `sample` statement that uses another prior, coupled with an `obs` argument. In the example below we first sample cutpoints `c_y` from the `ImproperUniform` distribution with `constraints.ordered_vector` as before, and then `sample` a dummy parameter from a `Normal` distribution while conditioning on `c_y` using `obs=c_y`. Effectively, we've created an improper / unnormalized prior that results from restricting the support of a `Normal` distribution to the ordered domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [00:03<00:00, 256.41it/s, 7 steps of size 5.31e-01. acc. prob=0.92]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "   b_X_eta      1.23      0.31      1.23      0.64      1.68    501.31      1.01\n",
      "    c_y[0]     -0.24      0.34     -0.23     -0.76      0.38    492.91      1.00\n",
      "    c_y[1]      1.77      0.40      1.76      1.11      2.42    628.46      1.00\n",
      "\n",
      "Number of divergences: 0\n"
     ]
    }
   ],
   "source": [
    "def model2(X, Y, nclasses=3):\n",
    "    b_X_eta = sample(\"b_X_eta\", Normal(0, 5))\n",
    "    c_y = sample(\n",
    "        \"c_y\",\n",
    "        ImproperUniform(\n",
    "            support=constraints.ordered_vector,\n",
    "            batch_shape=(),\n",
    "            event_shape=(nclasses - 1,),\n",
    "        ),\n",
    "    )\n",
    "    sample(\"c_y_smp\", Normal(0, 1), obs=c_y)\n",
    "    with numpyro.plate(\"obs\", X.shape[0]):\n",
    "        eta = X * b_X_eta\n",
    "        sample(\"Y\", OrderedLogistic(eta, c_y), obs=Y)\n",
    "\n",
    "\n",
    "kernel = NUTS(model2)\n",
    "mcmc = MCMC(kernel, num_warmup=250, num_samples=750)\n",
    "mcmc.run(mcmc_key, X, Y, nclasses)\n",
    "mcmc.print_summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proper Prior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If having a proper prior for those cutpoints `c_y` is desirable (e.g. to sample from that prior and get [prior predictive](https://en.wikipedia.org/wiki/Posterior_predictive_distribution#Prior_vs._posterior_predictive_distribution)), we can use [TransformedDistribution](http://num.pyro.ai/en/stable/distributions.html#transformeddistribution) with an [OrderedTransform](http://num.pyro.ai/en/stable/distributions.html#orderedtransform) transform as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [00:04<00:00, 244.78it/s, 7 steps of size 5.54e-01. acc. prob=0.93]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "   b_X_eta      1.40      0.34      1.41      0.86      1.98    300.30      1.03\n",
      "    c_y[0]     -0.03      0.35     -0.03     -0.57      0.54    395.98      1.00\n",
      "    c_y[1]      2.06      0.47      2.04      1.26      2.83    475.16      1.01\n",
      "\n",
      "Number of divergences: 0\n"
     ]
    }
   ],
   "source": [
    "def model3(X, Y, nclasses=3):\n",
    "    b_X_eta = sample(\"b_X_eta\", Normal(0, 5))\n",
    "    c_y = sample(\n",
    "        \"c_y\",\n",
    "        TransformedDistribution(\n",
    "            Normal(0, 1).expand([nclasses - 1]), transforms.OrderedTransform()\n",
    "        ),\n",
    "    )\n",
    "    with numpyro.plate(\"obs\", X.shape[0]):\n",
    "        eta = X * b_X_eta\n",
    "        sample(\"Y\", OrderedLogistic(eta, c_y), obs=Y)\n",
    "\n",
    "\n",
    "kernel = NUTS(model3)\n",
    "mcmc = MCMC(kernel, num_warmup=250, num_samples=750)\n",
    "mcmc.run(mcmc_key, X, Y, nclasses)\n",
    "mcmc.print_summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principled prior with Dirichlet Distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is non-trivial to apply our expertise over the cutpoints in latent space (even more so when we are having to provide a prior before applying the OrderedTransform). \n",
    "\n",
    "Natural inclination would be to apply Dirichlet prior model to the ordinal probabilities. We will follow proposal by M.Betancourt ([1], Section 2.2) and use [Dirichlet](http://num.pyro.ai/en/stable/distributions.html#dirichlet) prior model to induce cutpoints indirectly via [SimplexToOrderedTransform](http://num.pyro.ai/en/stable/distributions.html#simplextoorderedtransform). \n",
    "This approach should be advantageous when there is a need for strong prior knowledge to be added to our Ordinal model, eg, when one of the categories is missing in our dataset or when some categories are strongly separated (leading to non-identifiability of the cutpoints). Moreover, such parametrization allows us to sample our model and conduct prior predictive checks (unlike `model1` with `ImproperUniform`).\n",
    "\n",
    "We can sample cutpoints directly from `TransformedDistribution(Dirichlet(concentration),transforms.SimplexToOrderedTransform(anchor_point))`. However, if we use the Transform within the `reparam handler` context, we can capture not only the induced cutpoints, but also the sampled Ordinal probabilities implied by the `concentration` parameter. `anchor_point` is a nuisance parameter to improve identifiability of our transformation (for details please see [1], Section 2.2)\n",
    "\n",
    "Please note that we cannot compare latent cutpoints or b_X_eta separately across the various models as they are inherently linked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will apply a nudge towards equal probability for each category\n",
    "# (corresponds to equal logits of the true data generating process)\n",
    "concentration = np.ones((nclasses,)) * 10.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [00:05<00:00, 193.88it/s, 7 steps of size 7.00e-01. acc. prob=0.93]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                 mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "    b_X_eta      1.01      0.26      1.01      0.59      1.42    388.46      1.00\n",
      "     c_y[0]     -0.42      0.26     -0.42     -0.88     -0.05    491.73      1.00\n",
      "     c_y[1]      1.34      0.29      1.34      0.86      1.80    617.53      1.00\n",
      "c_y_base[0]      0.40      0.06      0.40      0.29      0.49    488.71      1.00\n",
      "c_y_base[1]      0.39      0.06      0.39      0.29      0.48    523.65      1.00\n",
      "c_y_base[2]      0.21      0.05      0.21      0.13      0.29    610.33      1.00\n",
      "\n",
      "Number of divergences: 0\n"
     ]
    }
   ],
   "source": [
    "def model4(X, Y, nclasses, concentration, anchor_point=0.0):\n",
    "    b_X_eta = sample(\"b_X_eta\", Normal(0, 5))\n",
    "\n",
    "    with handlers.reparam(config={\"c_y\": TransformReparam()}):\n",
    "        c_y = sample(\n",
    "            \"c_y\",\n",
    "            TransformedDistribution(\n",
    "                Dirichlet(concentration),\n",
    "                transforms.SimplexToOrderedTransform(anchor_point),\n",
    "            ),\n",
    "        )\n",
    "    with numpyro.plate(\"obs\", X.shape[0]):\n",
    "        eta = X * b_X_eta\n",
    "        sample(\"Y\", OrderedLogistic(eta, c_y), obs=Y)\n",
    "\n",
    "\n",
    "kernel = NUTS(model4)\n",
    "mcmc = MCMC(kernel, num_warmup=250, num_samples=750)\n",
    "mcmc.run(mcmc_key, X, Y, nclasses, concentration)\n",
    "# with exclude_deterministic=False, we will also show the ordinal probabilities sampled from Dirichlet (vis. `c_y_base`)\n",
    "mcmc.print_summary(exclude_deterministic=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
